program example_z_MySparseOperator
    use MySparseOperator
    use common_functions
    implicit none
	complex(kind=8), dimension(:), allocatable :: a
    integer, dimension(:), allocatable :: rowind, colind
    complex(kind=8), dimension(:), allocatable :: b, x

    type(mso_zsparse_matrix) :: Amtx, Bmtx, Cmtx

    integer :: msglvl
    integer :: maxiter(2), itercount
    real(kind=8) :: maxrelres, residual(2)
    integer :: isfull = 1

! create a sparse matrix by indices of elements
    write(*,*) '============= create a sparse matrix by indices of elements =============='
    if ( allocated(a) ) deallocate(a) 
    allocate(a(10))
    a = dcmplx( (/1:9,3/), (/-10:-1/) )
    call mso_sparse_create_coo( (/ 1:9, 4 /), (/ 1:9, 4 /), a, Amtx, 9)
    call mso_sparse_show(Amtx, '(f6.1)', isfull)

! add two sparse matrix
    write(*,*) '============= add two sparse matrix ======================================'
    if ( allocated(a) ) deallocate(a) 
    allocate(a(8))
    a = dcmplx(3.0d+0,5.0d+0)
    call mso_sparse_create_coo( (/ 1:8/), (/ 2:9/), a, Bmtx, 9)
    call mso_sparse_add(Amtx,dcmplx(1.d0),Bmtx,Cmtx)
    a = dcmplx(6.0d+0,7.0d+0)
    call mso_sparse_create_coo( (/ 1:7/), (/ 3:9/), a(1:7), Bmtx, 9)
    call mso_sparse_add(Cmtx,dcmplx(1.d0),Bmtx,Amtx)
    a = dcmplx(7.0d+0,-1.0d+0)
    call mso_sparse_create_coo( (/ 1:6/), (/ 4:9/), a(1:6), Bmtx, 9)
    call mso_sparse_add(Amtx,dcmplx(1.d0),Bmtx,Cmtx)
    call mso_sparse_show(Cmtx, '(f6.1)', isfull)
    write(*,*) '=========================================================================='
    call mso_sparse_add(Cmtx,dcmplx(1.d0),Cmtx,Amtx,'t')
    call mso_sparse_show(Amtx, '(f6.1)', isfull)

! delete two rows and two columns from a sparse matrix
    write(*,*) '====== delete two rows and two columns from a sparse matrix =============='
    write(*,*) '=== delete two rows ================='
    call mso_sparse_deleteRowCol(Amtx,(/2,3/),'R',Bmtx)
    call mso_sparse_show(Bmtx, '(f6.1)', isfull)
    write(*,*) '=== delete two columns ================='
    call mso_sparse_deleteRowCol(Bmtx,(/5,7, 10, 11/),'C',Cmtx)
    call mso_sparse_show(Cmtx, '(f6.1)', isfull)

! get the upper triangular part of a sparse matrix
    write(*,*) '========= get the upper triangular part of a sparse matrix ==============='
    write(*,*) '====== upper =  0 ==========='
    call mso_sparse_uppertriangle(0,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, '(f6.1)', isfull)
    write(*,*) '====== upper =  1 ==========='
    call mso_sparse_uppertriangle(1,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, '(f6.1)', isfull)
    write(*,*) '====== upper =  -1 ==========='
    call mso_sparse_uppertriangle(-1,Amtx,Bmtx)
    call mso_sparse_show(Bmtx, '(f6.1)', isfull)

! compute matrix - vector product of a sparse matrix
    write(*,*) '======= compute matrix - vector product of a sparse matrix ==============='
    if ( allocated(x) ) deallocate(x)
    allocate(x(6))
    x = dcmplx(2.d0,4.d0)
    call mso_sparse_deleteRowCol(Amtx,(/1,4,9/),'C',Bmtx)
    write(*,*) '======= the matrix is:'
    call mso_sparse_show(Bmtx, '(f6.1)',isfull)
    write(*,*) '======= the vector is:'
    call prtVec(x,'(f6.1)')
    call mso_sparse_mv(Bmtx,x,b)
    write(*,*) '======= the result is:'
    call prtVec(b,'(f6.1)')
    call cf_allocate(x,Bmtx%rows)
    x = dcmplx(2.d0,4.d0)
    write(*,*) '======= a new vector is:'
    call prtVec(x,'(f6.1)')
    call mso_sparse_mv(Bmtx,x,b,'c')
    write(*,*) '======= the result is:'
    call prtVec(b,'(f6.1)')


! use LU method (mso_dpardiso) to caculate the solution of a sparse linear system
    write(*,*) '======= caculate the solution of a sparse linear system ================='
    if ( allocated(x) ) deallocate(x)
    allocate(x(Amtx%rows))
    x = cmplx(3.3,-5.3)
    call mso_sparse_mv(Amtx,x,b)
    write(*,*) '================expected solution is:'
    call prtvec(x,'(f5.1)')
    call mso_sparse_pardiso(Amtx,b,x)
    write(*,*) '======= computed solution is:'
    call prtVec(x,'(f5.1)')

    call mso_sparse_mv(Amtx,x,b)
    call mso_sparse_setype(Amtx,3) ! set the matrix to be structurally symmetric
    call mso_sparse_pardiso(Amtx,b,x)
    write(*,*) '======= computed solution is:'
    call prtVec(x,'(f5.1)')

    call mso_sparse_add(Amtx,dcmplx(1.d0),Amtx,Bmtx,'C')
    call mso_sparse_mv(Bmtx,x,b)
    call mso_sparse_setype(Bmtx,-4) ! set the matrix to be Hermitian indefinite
    call mso_sparse_pardiso(Bmtx,b,x)
    write(*,*) '======= computed solution is:'
    call prtVec(x,'(f5.1)')










    if ( allocated(a) ) deallocate(a)
    if ( allocated(b) ) deallocate(b)
    if ( allocated(x) ) deallocate(x)
    call mso_sparse_destroy(Amtx)
    call mso_sparse_destroy(Bmtx)
    call mso_sparse_destroy(Cmtx)


end program example_z_MySparseOperator